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ABSTRACT 

We construct analytic models of incompressible, uniformly rotating stars in post-Newtonian 
(PN) gravity and evaluate their stability against nonaxisymmetric bar modes. We model the 
PN configurations by homogeneous triaxial ellipsoids and employ an energy variational principle 
to determine their equilibrium shape and stability. The spacetime metric is obtained by solving 
Einstein's equations of general relativity in 3+1 ADM form. We use an approximate subset of 
these equations well-suited to numerical integration in the case of strong field, three dimensional 
configurations in quasi-equilibrium. However, the adopted equations are exact at PN order, 
where they admit an analytic solution for homogeneous ellipsoids. We obtain this solution for 
the metric, as well as analytic functional for the conserved global quantities, M, Mq and J. 

We present sequences of axisymmetric, rotating equilibria of constant density and rest 
mass parametrized by their eccentricity. These configurations represent the PN generalization 
of Newtonian Maclaurin spheroids, which we compare to other PN and full relativistic 
incompressible equilibrium sequences constructed by previous investigators. We employ the 
variational principle to consider nonaxisymmetric ellipsoidal deformations of the configurations, 
holding the angular momentum constant and the rotation uniform. We locate the point along 
each sequence at which these Jacobi-like bar modes will be driven secularly unstable by the 
presence of a dissipative agent like viscosity. We find that the value of the eccentricity, as well 
as related ratios like 0^/(7rpo) and T/|VF| (= rotational kinetic energy / gravitational potential 
energy), defined invariantly, all increase at the onset of instability as the stars become more 
relativistic. Since higher degrees of rotation are required to trigger a viscosity-driven bar mode 
instability as the stars become more compact, the effect of general relativity is to weaken the 
instability, at least to PN order. This behavior is in stark contrast to that found recently for 
secular instability via nonaxisymmetric, Dedekind-like modes driven by gravitational radiation. 
These findings support the suggestion that in general relativity nonaxisymmetric modes driven 
unstable by viscosity no longer coincide with those driven unstable by gravitational radiation. 

Subject headings: gravitation — relativity — instabilities — stars: neutron — stars: rotation 

1. Introduction 

The identification of nonaxisymmetric modes of instability in rapidly rotating equilibrium configurations 
is a classic problem. Although a considerable amount of work has been done in Newtonian theory, where 
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a number of results are well established (see Chandrasekhar 1969a, hereafter Ch69 and section 2 for a 
summary and references), only a few investigations have been carried out so far in the context of general 
relativity. The first three dimensional (3D) perturbation computations in full general relativity to identify 
nonaxisymmetric instabilities driven by gravitational radiation have been carried out recently by Stergioulas 
& Friedman (1997, hereafter SF) and Stergioulas (1997). A earlier numerical investigation of the effects of 
relativity on the viscosity-driven bar mode instability was presented by Bonazzola, Frieben & Gourgoulhon 
(1996, hereafter BFG) and these results have been recently corroborated by a more detailed analysis 
(Bonazzola, Prieben & Gourgoulhon 1997). SF solved the coupled set of perturbed field equations, modeling 
uniformly rotating stars by a polytropic equation of state (EOS) and adopting the Friedman & Schutz 
(1995) criterion for the onset of the nonaxisymmetric instability to gravitational radiation dissipation. 
According to this criterion, a nonaxisymmetric mode becomes unstable when its frequency, as measured by 
an observer at infinity, vanishes. They find that relativistic models are unstable to nonaxisymmetric modes 
for significantly smaller degrees of rotation than for corresponding Newtonian models. The destabilizing 
effect of relativity is most striking in the case of the to = 2 mode, which can become unstable even for soft 
polytropes of index 7 > 'ycrit = 1-77 (n < 1.3), while the critical index in Newtonian theory is 'jcrit = 2.238 
(Jeans 1919, 1928; James 1964). This behaviour is in agreement with results of a semianalytical, PN 
analysis previously presented by Cutler (1991) and Cutler & Lindblom (1992, see also Lindblom 1995, 
Yoshida & Eriguchi 1997). In particular. Cutler (1991) derived the PN corrections to the pulsational modes 
of uniformly rotating stars. The resulting expressions were then used by Cutler & Lindblom (1992) to 
evaluate the critical value of the star angular velocity, ^crit, where the frequency of the mode passes through 
zero. This is the critical value at which, in absence of viscosity, these modes are unstable to the emission of 
gravitational radiation. By considering the 7 = 2 polytropc. they found that PN effects lowers by up to 10 
% the values of flcrit, concluding that PN effects tend to make the gravitational radiation instability more 
important. However, as noted by SF, the results of BFG, who investigated the effects of general relativity 
on the TO = 2, viscosity-driven "bar" mode instability, seem to suggest the opposite effect. BFG's method 
consists of perturbing a stationary, axisymmetric configuration, obtained from a two dimensional numerical 
simulation, and retaining only the dominant terms in the nonaxisymmetric relativistic perturbation 
equations. BFG found that the critical adiabatic index for the bar instability becomes higher than the value 
of James as the configuration becomes increasingly relativistic. This behaviour suggests that relativistic 
effects tend to stabilize the configurations. Accordingly, SF concluded that, in general relativity, the point 
of onset of the viscosity-driven and the gravitational radiation-driven to = 2 modes may no longer coincide 
as they do in Newtonian theory, and that the effect of relativity seems to be very different in the two cases. 
The main improvement presented in Bonazzola, Prieben & Gourgoulhon (1997) over their previous study 
consists in the fully three dimensional treatment of the shift vector. With this more detailed analysis, the 
stabilization of relativistic configurations is not only confirmed, but also strongly enhanced. Both SF's and 
BFG's findings are based on a numerical solution of perturbation equations and represent the only attempts 
to solve the relativistic problem in 3D to date. 

In this paper we reconsider the problem of relativistic rotating equilibria and bar mode instabilities 
from an analytic point of view. Specifically, we extend the earlier Newtonian treatment of Lai, Rasio & 

Shapiro (1993a, hereafter LRS), which is based on triaxial ellipsoid models of rotating stars and an energy 
variational principle, to post-Newtonian (PN) gravitation. We restrict our discussion to incompressible, 
rigidly rotating bodies and neglect any deviation from the ellipsoidal shape in the equilibrium configuration. 
The ellipsoidal approximation is exact for rotating incompressible stars in Newtonian theory, but it is only 
approximate for PN configurations. However, our formalism allows us to derive the analytic functionals 
for the main global parameters characterizing a rotating configuration (total mass-energy, rest mass and 
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angular momentum). By applying a energy variational principle to these functionals, it is possible to 
construct equilibrium sequences of constant rest mass and to locate instability points along the sequence. 

Our PN analysis is carried out in the framework of a 3+1 ADM splitting of the metric (Arnowitt, 
Deser and Misner 1962). In this respect, our solution sets up and should provide an important test-bed 
calculation for future numerical studies of 3D relativistic rotating configurations. Numerical relativity in 
3D is only in its infancy, and lacks a large body of known solutions which it can attempt to reproduce. 
Part of the complication is related to the fact that, in general relativity, 3D stellar systems are usually 
fully dynamical due to the generation of gravitational waves. As a result, it is typically necessary to 
solve the full set of Einstein's equations to determine the behavior of a 3D system. There are many 
asymmetric systems, however, in which the intrinsic dynamical timescale is much shorter than the timcscale 
for dissipation due to gravitational waves. Such is the case for any weak field, slow velocity system, or 
for a strong-field, high velocity system that is only slightly perturbed from stationary equilibrium. It is 
thus possible to discuss "quasi-equilibrium" configurations for such objects. Examples include rotating 
neutron stars (NSs) either with small compaction M/ R or small departure from axisymmetry, and binary 
neutron stars prior to reaching the innermost stable circular orbit (Baumgarte et al. 1997b). Following 
the discussion by Wilson & Matthews (1989, 1995) and Wilson (1990), Cook, Shapiro and Teukolsky 
1992 (hereafter CST92) have provided a simplified set of 3+1, ADM equations which are well-adapted 
for studying quasi equilibrium relativistic systems. For strong-field objects, these equations yield precise 
solutions to the initial value (constraint) equations, and approximate instantaneous snapshots of the object 
as it evolves, due to the emission of gravitational waves, provided it does so slowly. In this paper, we adopt 
the CST92 equations which are well suited to future numerical studies of the fully nonlinear equations for 
strong- field, quasi-equilibrium sources. We solve these equations in the PN approximation; at this order, 
the CST92 equations are exact. In this paper, we focus on the viscosity-driven, secular instability with 
respect to the formation of a barlike structure in a PN incompressible, rotating star. We find that PN 
solutions are unstable at rotation rates greater than in the Newtonian limit. Thus, in the PN treatment, 
nonaxisymmetric instabilities driven by viscosity set a less stringent limit on the maximum angular velocity 
than suggested by Newtonian theory. 

From the observational point of view, the issue of investigating the maximum spin velocity for a NS is 
manyfold. Rotating configurations subject to nonaxisymmetric instabilities could become important sources 
of gravitational waves emission and represent possible candidates for detection by laser interferometers 
now under construction, like LIGO, GEO and VIRGO (see e.g. Bonazzola & Marck 1994; Thorne 1987; 
Lai & Shapiro 1995; Schutz 1997 and reference therein). There are a number of plausible astrophysical 
situations in which this effect can set in. Numerical simulations of coalescing binary neutron stars show 
that merger into a single object is the probable end point of close binary evolution (see e.g. Rasio & 
Shapiro 1992, 1994; LRS; Lai, Rasio & Shapiro 1993b, 1994a,b,c). The newly born object may be rotating 
quite rapidly (Baumgarte et al. 1997b) and might be driven secularly unstable. Core collapse in massive, 
evolved stars or accretion-induced collapse of white dwarfs, also has been suggested as observable sources 
of gravitational radiation (Lai & Shapiro 1995 and references therein). Coalescing white dwarf binaries are 
though to be progenitors of type la supernovae (Iben & Tutukov 1984; Yungelson et al. 1994), or, in certain 
cases, of isolated millisecond pulsars (Chen & Leonard 1993). In all these scenarios, the newly born NS can 
\mdergo a secular or dynamical evolution, breaking its axisymmetry, provided that its spin velocity is larger 
than the critical value. Alternatively, accretion onto NSs in X-ray binary systems or in Thorne-Zytkow 
objects can spin up the compact star until reaching the critical value of T/|VF|, at which point the 
Chandrasekhar-Priedman-Schutz (CPS) instability (Chandrasekhar 1970; Friedman & Schutz 1978) may 
drive a nonaxisymmetric mode unstable, powering radiation. This "forced gravitational emission", is 
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particiilaiiy appealing for maintaining steady, periodic emission, with the total amount of accreted angular 
momentum balancing the amount that is radiated away via gravitational radiation (Wagoner 1984; Schutz 
1997). 

We point out that our analytical results assume incompressible and rigidly rotating bodies, and both 
assumptions have been introduced mainly to make an analytical treatment tractable, to PN order. However, 
for applications to NSs, we note that both assumptions are not entirc^ly "ad hoc". Viscosity tends to drive 
NS to rigid rotation (Fricdmann & Ipser 1987). In addition, recent many body calculations suggest that 
the EOS of dense nuclear matter is relatively stiff (see, e.g. Wiringa, Fiks & Fabronici 1988 and references 
therein). In view of these considerations, application of our results to realistic situations may not be 
unreasonable, at least as a first approximation. 

In section 2, we review the physical problem and summarize the results of previous investigations. 
The variational principle of LRS is reviewed in section 3, while in section 4 wc set up the mathematical 
problem in the PN approximation. We assemble the relevant set of 3+1 equations for the gravitational field 
and present analytic solutions for both the metric coefficients and the global conserved quantities. These 
results are then used to build equilibrium sequences in section 5 and locate the bar mode instability point 
in section 6. Discussion and conclusions follow in section 7. 



2. Previous Investigations 

The problem of the equilibrium shape and stability of an incompressible, rigidly rotating configuration 
admits an exact analytic solution in Newtonian gravitation (Ch69) . In this case the rotating configuration 
takes in axisymmetry the equilibrium shape of a Maclaurin spheroid. However, nonaxisymmetric instabilities 
can develop in rapidly spinning spheroids when the ratio T/|VF| of the rotational kinetic to the gravitational 
potential energy becomes sufficiently large. At the critical value = 0.1375 the equilibrium sequence 

of MacLaurin configurations bifurcates into two other branches of triaxial equilibria, the Jacobi and the 
Dedekind ellipsoids. Since the Maclaurin spheroids are dynamically unstable only for > 0.2738, 

the bifurcation point is dynamically stable. However, in presence of a suitable dissipative mechanism such 
as viscosity or gravitational radiation, this point becomes secularly unstable to a / = 2 m = 2 bar mode. 
Secular instability may play an important role in limiting the maximum rotation of neutron stars. Viscosity 
dissipates the mechanical energy E, but preserves angular momentum J. Consequently, a Maclaurin 
spheroid undergoing a viscosity driven instability terminates its evolution as a Jacobi ellipsoid, with a lower 
value of E but the same value of J (and rest mass Mq). The Jacobi solution is rigidly rotating, so that the 
viscous dissipation stops once it is formed. Alternatively, the Dedekind ellipsoid represents a lower energy 
state (with respect to the Maclaurin solution) for a given circulation. A rapidly spinning configuration can 
be unstable to the emission of gravitational radiation at the bifurcation point (the CFS instability). Since 
this process does not conserve angular momentum but conserves circulation, under a CFS instability the 
growth of the bar mode leads to the deformation of a Maclaurin spheroid into a Dedekind ellipsoid. At this 
final state, which is stationary, the emission of gravitational waves stops. The competition between the 
Jacobi-like and the Dedekind-like modes is governed by the ratio of the strength of the viscous stress to the 
gravitational radiation reaction force, and this ratio depends crucially on the internal properties of the star. 
Although the determination of a realistic equation of state for NSs and, as a consequence, of the strength of 
the viscosity, is still an open issue, viscosity is thought to be more important for old NSs and gravitational 
radiation in the hottest, newly born objects (see e.g. BFG and references therein). The situation is more 
complicated when both viscosity and gravitational radiation act together, since they tend to cancel each 
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other, stabilizing the star (Lindblom & Detweiler 1977; Lai & Shapiro 1995). 

A number of efforts in Newtonian physics have been devoted to the extension of the Ch69's results 
to more realistic, compressible fluids, modeled by a polytropic equation of state. For rigidly rotating 
polytropes, bifurcation to triaxial configurations can only exist when the adiabatic index exceeds a critical 
value, 'fcrit = 2.238 (Jeans 1919, 1928; James 1964). This is because the EOS must be stiff enough to make 
the angular velocity at the bifurcation point lower than the limiting value Clk at which the centrifugal 
force balance the gravitational force at the equator (mass shedding limit). Moreover, Ipser & Managan 
(1985) demonstrated that the m = 2 Jacobi-like bifurcation point and the m = 2 Dedekind-like point have 
the same location along uniformly rotating, polytropic sequences, as in the incompressible case (see also 
LRS; Lai & Shapiro 1995). Typically, rotating equilibrium stellar models must be constructed mmierically. 
Models have been constructed by a number of authors, using both polytropic (see e.g. Bodenheimer & 
Ostriker 1973; Ipser & Managan 1981; Hachisu & Eriguchi 1982; Hachisu 1986a,b; see Tassoul 1978 for an 
extensive set of references) and more realistic equations of state for both white dwarfs and NSs (see e.g. 
Ostriker & Tassoul 1969; Durisen, 1975; Hachisu 1986a; BEG and references therein). 

LRS constructed triaxial ellipsoid models of rotating polytropic stars in Newtonian gravity, using an 
ellipsoidal energy variational method. This approach was originally introduced by Zel'dovich & Novikov 
(1971, see also Shapiro & Teukolsky 1983, hereafter ST) in the axisymmetric case, to investigate the 
stability of a polytropic star against gravitational collapse. The main advantage of the method comes from 
its simplicity: all results are analytic or quasi-analytic, lending themselves to straightforward physical 
interpretation. This is in part because the method deals directly with global, conserved quantities. When 
quantities like the total energy and the total angular momentum are determined along an equilibrium 
sequence, the evolution of the system can be tracked as it loses E or ,/ by some quasi-static dissipative 
process. Generalizing the same approach to the triaxial case, LRS were able to construct equilibrium 
sequences for compressible analogues of most classical incompressible objects, like isolated Maclaurin, 
Jacobi, Dedekind and and Riemann ellipsoids and binary Roche, Darwin and Roche-Riemann ellipsoids. 

Neutron stars are relativistic objects and the analysis of their equilibrium and stability must be 
based necessarily on general relativistic models. The structure of a rotating axisymmetric star in general 
relativity has been investigated numerically by a number of authors (see e.g. Butterworth & Ipser 1976, 
hereafter BI; Friedman, Ipser & Parker 1986; CST92; Cook, Shapiro & Teukolsky 1994a,b; Cook, Shapiro 
& Teukolsky 1996, hereafter CST96 and references therein), but the first fully relativistic computations 
of nonaxisymmetric instabilities have been presented only very recently (SF; Stergioulas 1997). Erom the 
analytical point of view, even less is known. An exact treatment of the radial oscillations of a gaseous mass 
in general relativity is possible (Chandrasekhar 1964), but a similar analysis of nonradial oscillations is not 
to be expected. In fact, apart from the difficulties associated with the solution of the Einstein equations 
without any presupposed symmetry, allowance must be made in the relativistic regime for the emission of 
gravitational waves. However, some insight into the nature of the general relativistic effects can be obtained 
by examining the problem in the PN approximation, i.e. at a level in which gravitational radiation plays 
no role. PN effects on the equilibrium of uniformly rotating, homogeneous bodies have been extensively 
investigated in a series of paper by Chandrasekhar (1965a, b, 1967a, b,c, 1969b: see also Chandrasekhar & 
Nutku 1969, hereafter CN, for the PPN corrections) using the tensor virial formalism. Whenever possible 
in this paper, we make a direct comparison between our expressions and the corresponding PN results 
derived by Chandrasekhar. In particular, starting from the PN equations of hydrodynamics, he derived 
the equilibrium relation between the eccentricity and the angular velocity. Although he obtained integral 
expressions for the global conserved quantities, he did not evaluate them or give explicit formulae for the 



- 6 - 



PN corrections to the rest-mass, angular momentum and binding energy. A different method, based on the 
solution of the PN Poisson equation in oblate spheroidal coordinates, was presented by Bardeen (1971). 
The formalism turns out to be simpler with respect to the Chandrasekhar one, but this approach does not 
allow an immediate generalization to the triaxial case. In 3D, the general eigenfunctions resulting from the 
separation (Lame functions) arc indeed available, but many details about their general properties are not 
well studied. No analytic investigations of the location of the secular stability point in general relativity 
were provided in these earlier studies or, to our knowledge, elsewhere in the literature. 



3. The Energy Variational Method 

In this section we briefly review the energy variational approach in Newtonian theory (LRS), to 
introduce the basic concepts that we use in our PN calculations. 

Consider a self gravitating, isolated system. Each configuration (whether or not in equilibrium) can be 
specified by the total energy E and a number of conserved global quantities, such as the rest mass Mq and 
the total angular momentum J. Since E can be always written in terms of the fluid density and velocity 
flelds p(x), v{x) 

i? = i?[p(x),z;(x);Mo,J,...] , (1) 

the equilibrium state can be determined by extremizing this functional with respect to all variations in 
both /9(x) and i'(x), under the constraint that the conserved quantities are unchanged. Direct application 
of such a variational approach to a multidimensional system is a computationally challenge task. However, 
as discussed by LRS, a great simplification arises when we can replace the infinite number of degrees of 
freedom containc^d in p{x) and v(x) by a limited number of free parameters ai,a2, ■ ■ ■■ This we can often do 
for sufflciently simple systems under suitable simplifying assumptions. The total mass energy then becomes 

E = E[ai,a2,...;Mo,J,...] (2) 

and the equilibrium configuration is determined by extremizing this functional according to 

dE 

— = i = 1, 2 . . . (equilibrium) (3) 

under the constraint that Mq, J, . . . arc conserved. 

The onset of the instability can be then determined from 

( o„ ) = *)i = l)2,..., (onset of instability) (4) 

9 



daidoj J p 



where the subscript ^eq' indicates quantities evaluated along the equilibrium sequence. Clearly, whether the 
instability actually arises depends on the presence of a suitable dissipative mechanism which preserves the 
conservation laws assumed in the construction of the equilibrium model. 

Now consider a simple application: a homogeneous, uniformly rotating, Newtonian fluid system, with 
density po and angular velocity fl. In the incompressible case, the internal energy vanishes and the total 
energy is given by 

E = T + W, (5) 

where T and W are the rotational kinetic and the gravitational contributions, respectively. Assume that 
the surface of the conflguration is ellipsoidal in shape. Then the geometry of the system is completely 
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specified by the values (j = 1 ... 3) of the three semiaxes of the outer surface, where the pressure vanishes. 
However, foUowing LRS, it is more convenient to introduce an equivalent set of parameters defined by 

/ \ 2/3 /X 2/3 

X.-^A ,A. = r^) , (6) 



ai / \a2 



and 



R = (aiflsas)'/' ■ (7) 



Note that R represents the radius of the spherical configuration with the same volume as the rotating one; 
it is not, in general, related to the equilibrium state. The gravitational potential and the kinetic energies 
can be written as 

3 GAfg Ich ^ SGMj 

5 R 2{a,a2asf' 5 i? ^ ' ^ ' 

72 t2 

T^ — ^—h, (9) 

where Mq — AirpoR'^ /3, J = ftl and Is = 2MqR^ /h is the momentum of inertia of a sphere of the same 
volume. The momentum of inertia of the ellipsoid is / = Is/h, and the two dimensionless ratios / and h 
are defined as 

2i?2 2\ \l A2 



/=777i = o(^ + ^+^3AlA2) , (10) 



2R^ 2\l\l ^^^^ 



Here Ich = Tli-^i^l [called / in Ch69, equation (3.15)], and the dimensionless coefficients Ai can be 
calculated in terms of standard incomplete elliptic integrals involving only the axis ratios [see equations 
(3.33)-(3.35) in Ch69]. In the spherical limit ai — a2 — — R and f = h = 1. 

The total mass-energy may now be written as 

3 A/f2 

E = ---^{IcH + th) , (12) 

where 

ilsMl Airpo h ' ^ ' 

The equilibrium sequence can be constructed by varying E according to 

while holding Mq, J constant. For an incompressible fluid po = constant and there no variations with 
respect to R. This gives 

= (1^2). 

The derivatives of h and Ich with respect to Ai and A2 are reported in LRS [equations (A3), (A9)]. 
Exploiting these formulas and using expression (p^), conditions (|l^) can be cast in the form 

= ^-alA,-^-Ic,-\^j2al~al) (16) 
= (1^2), 
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and these can be combined by adding each one to 1/2 times the other, yielding 



q2 ^ q2 



From Ch69, we write 



A, 



A 



B. 



B 



A., 



2 , A 2 



ofAij 



3Auaf + A.jaf + Aikaf {i ^ j ^ k) 



and add to each side of (17) the quantity 2a\a\Ai2^ obtaining 



2Si2 = a: 



2Si2 ) = 2 {a\alAi2 - a^^s) • 



The latter equalities allow a solution with ai 7^ 02 if and only if 



2 2a 2 a 

0^02^12 = 03^3 



2B, 



(Jacobi ellipsoid) 



(17) 



(18) 

(19) 
(20) 



(21) 

(22) 
(23) 



and, as it is well known, these two conditions determine the equilibrium sequence of Jacobi ellipsoids. In 
particular, the first represents a relation between the two axial ratios, while the latter gives the angular 
velocity. When ai = 02, we obtain 



q2 

All = 03^3 , = 2Sii . 



or 



^2 



— = 2 ( Ai - ^^3 



(Maclaurin spheroid) 



(24) 



(25) 



which recovers the Maclaurin sequence. To determine the condition for the onset of the secular instability 
to nonaxisymmetric perturbations (LRS), we solve 



det 



Since we have 



dXl 



= 0, t,j = l,2. 



d^E 



the determinant vanishes at the two points 

( d'^E 
\dXl 



± 



eq 



( d^E 
\dX1dX2 



(26) 



(27) 



(28) 



eq 



However, only the plus sign is relevant: it is easy to demonstrate that the minus sign corresponds to a 
negative value of the ratio /{ttpq) and must be discarded. The relevant solution gives the condition 



dxl 



t 



+ t 



8X18X2 



(29) 



- 9 - 



Finally, using the equations of the second derivatives reported in LRS [equations A(3), (A9)], equation 
can be cast in the form 



(H) 




(30) 



which is the exact expression for the secular instability point in Maclaurin spheroids with respect to 
nonaxisymmetric perturbations (see Ch69, LRS). Given the constraint that the deformations be ellipsoidal, 
we have identified the onset of the bar mode (to = 2) instability point. Here the instability we have 
located would be triggered by the presence of viscosity, since uniform rotation was assumed in building the 
equilibrium sequence and maintained in the variations, while holding fixed Mq and J (but not circulation). 
We can also recover the result that this condition occurs at the precise point where the Jacobi sequence 
bifurcates from the Maclaurin sequence [see (^3|) , (|4|)] . 

Using the variational method in the framework in Newtonian physics, LRS were able to investigate 
a number of problems. In particular, they constructed approximate hydrostatic equilibrium solutions for 
rotating polytropes, either isolated or in binary systems, and presented the compressible generalization 
of the classical sequences of Maclaurin spheroids, Jacobi, Dedekind and Riemann ellipsoids and Roche, 
Darwin and Roche-Riemann binaries. For the case of incompressible Maclaurin sequences the ellipsoidal 
approximation is exact and the results derived from the variational principle are exact (Ch69). In this paper 
we will use the same variational approach to investigate the ellipsoidal configurations in FN gravitation. 



We will construct a variational expression for the total mass-energy M of a rotating ellipsoid in 
FN gravitation. The variation in M is equivalent to a variation in i? = M — Mq for fixed Mq, and we 
will consider incompressible configurations with constant rest mass density pQ. For a FN treatment the 
evaluation of all the integral quantities (e.g. M, Mq and J) requires a knowledge of the metric, which 
must be determined self-consistently with the matter profile by solving the Einstein field equations. The 
field equations form a set of coupled, nonlinear partial differential equations and, in general, must be 
solved numerically. The solution of these equations for a rotating star in general relativity has been 
tackled by many authors (see e.g. CST92; Cook, Shapiro & Teukolsky 1994a, b; CST96, and references 
therein). A great simplification arises if we adopt the "conformal approximation" of CST96 (see also 
Wilson & Matthews 1989, 1995; Wilson 1990, Wilson, Matthews & Marronetti 1996). This form for the 
metric is only approximate (although the approximation is very good) for rotating stars in full general 
relativity (CST96), for two main reasons. First, even in absence of gravitational radiation (e.g., the case 
for stationary, axisymmetric equilibrium rotating stars) the exact solution is not conformal. However, as 
shown by CST96, the deviation from conformal flatness is small (< 1% even for highly relativistic stars). 
The conformal approximation is exact in FN gravitation and departures only enter at the 2 FN order. 
The second reason is that emission of gravitational radiation causes the system to evolve in time, so the 
conformal decomposition of CST96 does not hold in general. However, the timescale for the evolution due 
to gravitational radiation is much longer than the dynamical timescale, and in many situations relativistic 
systems can be treated to be in quasi-equilibrium (see e.g. Baumgarte et al. 1997a, b). An analogous 
approximation is often used in stellar evolution calculations, where the relevant evolution timescales are 
the nuclear or Kelvin-Hclmholtz timescales, so that the stars maintain (quasi) hydrostatic equilibrium on 
a dynamical timescale. Gravitational radiation only enters at the 2.5 FN order. Henceforth, we restrict 
our attention to configurations in FN theory, where the 3-metric is conformally flat. We choose to use the 
3+1 decomposition of Einstein's equations in our analysis, rather than the conventional Chandrasekhar 



4. 



The Post Newtonian Solution 
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(1965a) PN expansion. We do so in preparation for future numerical treatments of nonaxisymmetric 
instabilities and their nonlinear evolution, which are best tracked in 3+1 form (see also BonazzoUa, Frieben 
and Gourgoulhon 1997). 

Henceforth we confine our attention to the case of incompressible, rigidly rotating models, for which 
the PN system can be solved analytically provided we adopt a ellipsoidal model for the matter profile. 

We solve the conformal ADM equations for the metric in this section. We introduce our approximations, 
section 4.1, referring to CST92 and CST96 for the notation and for a more detailed discussion. The final 
metric is presented in section (4.2) and expressions for the conserved quantities are derived in (4.3). Our 
field equations are equivalent to those of CN in the PN limit. 



4.1. Basic Equations 

Let us consider an isolated, self-gravitating, homogeneous system with rest mass density po and total 
mass-energy density p. Assume that it is uniformly rotating with a constant angular velocity Q. In the 
incompressible limit, the internal energy is zero, and the energy-momentum tensor takes the form 

T^. = {po + P) U^U, + Pg^, , (31) 

where t/^ is the fluid four-velocity, i?^^ are the metric coefficients and P is the pressure. Here Greek indices 
p,v, . . . range over ... 3, while Latin indices i, j, . . . range over 1 ... 3; geometrized units {c — G = 1) are 
used throughout. We assume that the outer surface (where p = 0) is a triaxial ellipsoid with scmiaxes ai, 
«2, «3-0 Following CST96, we start from the most general expression for the metric in a 3+1 form 

ds^ = -a'^dt^ + {dx' + /3'dt) [dx^ + pHt) , (32) 

where a and /3* are the lapse and the shift functions, respectively. We choose a conformally flat 
decomposition of the spatial metric 

= • (33) 

Here is the conformal factor and fij is the Euclidean metric in the adopted coordinate system. In the 
following we will use cartesian coordinates Xi, i — I . . .3. Following CST96, we set dt (7~^^^7ij) — and 
adopt the maximal slicing condition K — Kl = 0, where Kij is the extrinsic curvature. We obtain [see 
equation (4) by CST96] 

K^, = ^ (d,P, + D,p, - ^7.,-Dfc/3' ) • (34) 

Here Di indicates the covariant derivative with respect to 7ij. The metric coefficients depend on the 
three functions a, /9* and 5" that, following CST96, can be derived as a solution of a system of ADM 
partial differential equations [see equations (8), (14) and (18) in CST96]. The first two equations are the 
Hamiltonian constraint equation and the lapse equation {dtK — 0) 

V^^- = 27r*'V, (35) 



^Here the ellipsoid is defined in coordinate space, so that the proper shape of the PN configuration is not in general ellipsoidal. 
Although convenient mathematically, such a trial function is only an approximation (except in the Newtonian limit, where it is 
exact) . 
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l^^K.jK'^ + 2Tr-^^ [p + 2S) 



(36) 



while the differential equation for the shift vector can be obtained by substituting expression ( |34| ) into the 
momentum constraint 

D,K'' = SvrS*^' . (37) 

The final result is 

V2/J^ + iv* (V,/3^) = [U/,a ~ |v,*^ (^V^73^ + V'/3^" - ^r-' Vfe/3'=^ (38) 

In the previous expressions and denote the flat space covariant derivative and the Laplacian operator. 
The density p appearing into equations (^) , ( ^ ) can be derived from the stress-energy tensor 

p = n'^n'^T,^ - (po + P) (aC/*) ' - P , (39) 

where is the normal vector to a f = constant surface. Note that p = po for a nonrotating sphere, since 
there is no internal energy. However, the two densities differ in general and p is not a constant, due to 
rotational energy contributions. The source term S and the momentum source 5* are given by 

S = Y^T.j = (po + P) [(«[/*)' - l] + 3P , (40) 

= -7>fer^'' = (po + P) {aU') f'U, . (41) 

Equation ( p8[ ) can be conveniently reduced to two simpler equations by introducing the decomposition 

= G' ~ iv'P . (42) 

The two equations that must be solved now become 



V^G^ = ( iv^a - ■^Vj^' ) ( + V'(3^ - '^P^kP'" ) + 167ra*'*5' , (43) 



V^P VfeC-' . (44) 



In the fully relativistic case the solution of the coupled ADM equations (^5|), (|36|), ( ^431) and ( [44|) 
represents a nontrivial problem and must be tackled numerically. In this paper we work at the PN order 
and, in the ellipsoidal approximation, the equations can be solved analytically. The metric and the stress 
tensor will be expanded as sums of terms of successively higher order in the expansion parameter 1/c^, 
while each ADM equation will be decomposed into a series of equations of successively highly order in 1/c^. 
The first PN correction will refer to the terms that are 0(c~^) (i.e., 0{M/R), where i? is a length scale of 
the problem) higher than the corresponding Newtonian terms in this expansion. We do not restrict our 

analysis to slow rotation (Hartle 1967), whereby one requires \J^VL <^ 1. In that context, rotation is 
considered "slow" if its effects on the structure of the star are relatively small. Here we allow arbitrary fast 
rotation, so that VP' is permitted to reach ^ {M/ R^^ and stars can suffer considerable rotational distortion. 



To obtain equations correct to the PN order, in the right-hand sides of (|35|), ( |36| ) and ( p8D we need 
retain only the contributions of order 

PO, Po-^- (45) 
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Consider the Hamiltonian constraint and the lapse equations. From equations (|3^) and (|4l|) we have, at the 
leading order, 

S' ~ pov ^ - S'R^ - povR^ - , (46) 

R 



which yields [see equation (| 



2 , ^2 



K.K'^-Vld -Po^- (47) 



This means that, at our order of approximation, we can safely drop all terms involving the extrinsic 
curvature from equations ( p5[ ) and (|36|), which then reduce to the simpler form 

V^* = -27r*V, (48) 

(a^-) = 27ra*5 (p + 25*) . (49) 

We thus need approximate expressions for the source terms p and S, and these can be obtained by 
expanding the product {aU*) appearing in ( |39| ) and (|4^). The normalization condition UuU'^ — —1 yields 

{aU'f = l^fW^Uj. (50) 

We introduce = and consider a velocity field corresponding to uniform rotation with angular 

velocity Vl = [/"^/C/* about the X3 direction. Then 

Vi = -Q.X2 , V2 = ^xi , U3 = , (51) 

and — Vl^{x\ + x"^). At the leading term, it follows that 

(aC/*)'«l + w2 + ... . (52) 

Substituting back in the definition of p and S and recalling that P/po ^ 0{M/R), we find that, at 
Newtonian order, these two quantities can be approximated as 

p^Po{l + v^ + ...), (53) 

S^Po(^v^ + 3^ + ..}j . (54) 

We also need the pressure distribution in the equilibrium configuration. Terms involving P appear only in 
the PN correction. This means that we can use the Newtonian result for P (see e.g. Chandrasekhar 1965a, 
hereafter Ch65a) 



P 

— = T^PO 
PO 



l^' + A) ■ (55) 



The functions 4' and a now can be derived by solving the two equations ( ^ ) and (^). The conformal 
factor can be expanded as: 

4- = 1 - $/2 , (56) 

where $ = + $ pat and where 

V^^AT EE 47rpo- (57) 

Note that <i>jv represents the Newtonian potential and coincides with the quantity in the Chandrasekhar's 
notation. Linearizing 5''^ sa (1 — 5<f>/2) and substituting (|53|), (p6|), (57) into (^8|), we obtain 

V^$PAr = -lOTrpo^'Af + 47rpot^^ • (58) 
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2 


2 



In a similar way, the lapse equation can be linearized by introducing the expansion = 1 + 0, where 
= Sat + 6pAr and V^Oat = 27rpo, so that Sat = <f>Ar/2. Writing 

*"«(l-^...)'«l-2$Ar..., (59) 

and 

p + 25* w po (^1 + + 6^ . . , (60) 

we obtain 

V^GpAT = 6^po K + 2— - ]:^n] ■ (61) 
The corresponding expansion of the lapse function is 

a = (i + e)vi/-i«i + ^ + ^ + ^ + eA. + eA.^ + ep^. (62) 

It is straightforward to identify the Newtonian and the PN contributions to a as 

(63) 

apN = ^ + ^ + QpN ■ (64) 

Parenthetically, we note that a at and apjv coincide with the correspondent terms derived by CN (see 
Appendix A. 2). Summarizing, to derive the PN corrections to the lapse function and the conformal factor, 
we need to solve the system of elliptic equations (|57|), ( |58| ) and (|6l|). 

Let us now focus on the shift vector, and introduce the expansions 

=/5^^ + 0(c-5), ^ G'pj, + 0{c-'') , B^BpN + 0{c-'>). (65) 

At our order of approximation, we only need retain the leading terms in these expansions, of order 
~ vM/R ^ c^^. As a consequence, we can legitimately drop the first term on the right hand of equation 
(^3|), being of order 

M 13' 

B^R^ (66) 

(note that Va - V$ M/R^). This yields 

V^G'pN w lQTTa<^^S' « 167r5* , (67) 
and we are left to derive the leading term in 5*. In expression ( pl| ) we have 

l^'U, = i^'gjuU'' (68) 

which yields to leading order 

« p^v' . (69) 
By substituting expressions ( ^l| ) for the components of v', we can finally write the two equations 

V^G^^ = 167^po^'^ = -167rpof^a;2 , (70) 
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V^G^^ = IGttpow^ = Wnponxi , (71) 
while GpjY = 0. Once Gpjy is known, Bp^ can be obtained simply by solving 

Y72n dGp^ 9Gpj^ 

and /3pjY will follow from (p^). 

As it will be shown in the next section, the full system of elliptic equations (|5^), (|6l|), (|7C|)-(|7^) 
admit an analytical solution (up to well known elliptic integrals). 



4.2. T/ie Metric: Analytic Solution 

The Newtonian potential at an internal point Xi of a homogeneous, ellipsoidal configuration with 
semi-axes has the well known analytical form 



= -Po / 1 — ^rf^x' = -Trpo I ^c/i - > ' A.x^ I , (73) 



[see e.g. Ch69, equation (3.40)]. This function provides the leading term in the expansions for and a 
[equations (^6|), (|63|)]. Hence we only need to determine the post-Newtonian corrections ^pn, otpN and 
/^PAT to obtain the metric at the PN order. The main advantage of our decomposition comes from the 
fact that the PN contributions can all be written in integral form by exploiting the Green's functions 
of the corresponding differential equations. Moreover, in our homogeneous, rigidly rotating case, the 
aforementioned integrals involve familiar quadratures of the kind 



The first one is simply — <i>Ar, while the others coincide with the two Newtonian potentials Di and Dij given 
by Ch69 [see equations (3.120), (3.131)]. As a result, we are able to derive explicit expressions for the three 
PN corrections ^pn, ctPN and /3p^ in terms of elliptic integrals. 

Consider first the PN correction to the conformal factor. Substituting = (.x^ + Xj) into the right 
hand of equation (^8|), the formal solution for ^pn can be expressed as 

2 J |x-x'| J |x-x'| 

Inserting the analytical expression ( [73|) for (f>jv in the integrand of ([75|), it is easy to recognize that ^pn 
can be written as 



5 . , 5 



= ^TTPalch'^N + ^vrpo V A,Du - n'^ {Dn + D22) , (76) 



where the potential Da, written in terms of the index symbols A^jj, and B^j^ , is [see equation (3.132) in 
Ch69, Appendix D] 



A» = P^ I TZT^f^' (77) 

A^^ - A,,jx] xi + ^npoai - 2 ^ By + ^ B.jux^^xl 
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The derivation of ap n and /3f ^ is carried out in a similar way. The formal solution of equation ( |6l| ) is 

and, by substituting <i>Ar and the pressure profile (|55|), we obtain 

GpAT = 37rpo (Asal + $Ar + ^rrpo ^^^^^ " 3^^' (^n + ^22) ■ (79) 

Once OpAT and ^pn are known, the PN correction ap^ follows trivially from (|6^). 

Consider next the shift vector. Solving equations (|70|), for Gpjy ^pw yields 

hence 

G},jv = 4r!D2 , G|jv = -417^1 , (81) 
where [see equation (3.121) by Ch69] 

A - Po / ^^rf'x' (82) 

Then, by differentiating expressions ( |8^ ) with respect to Xi, we obtain the source term appearing in the 
differential equation (|7^) for -Bpjv- This gives 

V^BpN — 87rripo^i2 (a? ~ 02) 2^1X2 , (83) 
where we used the symmetry property yli2 = A21. The solution is 

BpN = -2npoAi2 (al - al) J ^^^^'x' (84) 
= -2nAi2 {aj - aj) D12 , 

where 

XiX2- (85) 



D12 = TTpoalaj ^Ai2 - Ai2tx'^^ 



Note that in axisymmetry, when ai = a2, B ~ and the shift vector is divergence-free [see equations (p^), 
(^)]. To obtain the shift from (|4|), we must evaluate the gradient of the £'12 according to 

\/Wi2 = {5ii + Si2) - 2npoalalAi2iXiX2Xi , (86) 

Xi 

where 5ij is the Kronecker delta. The final result is 

1 - 

/3),^ = 40^2 + 7:^^^12^ ^ (i^i2 - 2^poa?a22;?a:2Ai2i) (87) 

2 xi 

1 - 

= -AnDi + -nAi2— ^{Di2-2npoalalxixlAi22) 

2 X2 

PpN = -^Ai2 (al - al) TTpQalalxiX2X3Ai23 . 
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Equations (^3|), (Jt^), ( [79| ) and ( ^7| ) completely determine the metric functions to PN order in the ellipsoidal 
approximation.^ We compare this solution in Appendix B.l to the exact solution for the interior metric 
in full general relativity for the case of a spherical, homogeneous configuration. We verify that our results 
agree with the PN expansion of the exact solution. 

We have checked our results whenever possible with the corresponding PN expressions derived by CN, 
and they agree. A comparison is provided in Appendix A. 



4.3. The Evaluation of the Conserved Quantities M , Mq and J 

4-3.1. The Integral Forms 

With the metric coefficients derived in section (4.2) we can now derive explicit expressions for the total 
mass-energy M, the total rest-mass Mq and the angular momentum J of a homogeneous ellipsoid. These 
functionals will be then employed to determine the properties of the equilibrium configurations, using the 
energy variational method. 

First consider the total baryon rest mass Mq, [CST92, equation (50)] 

Mo=Po f U*V^d^x, (88) 
Jv 

where g — det g^^ and V is the volume of the ellipsoid. In our adopted conformal gauge 

Mo = po / (aC/*) 'i'^d^x . (89) 



V 



In order to evaluate the PN correction to Mq, we need to know the integrand in (|89|) up to order M^/R^. 
To this order the quantity can be approximated as 



= 1^1- -j (90) 

w 1-3$ + — 
4 

« l-3$Ar-3$PAr + ^$?^ + 0(AfV^^) 

In addition, we need a similar expansion for the quantity at/* and this can be derived from the normalization 
condition (^0|). After some algebra, we obtain 

{aU*Y « i + (u*)^q,^[y^ + 2{pl^v2+PUvi)+---] (91) 

« 1 + ([/*)' (1 - 2$Ar + . . .) [v^ + 2 {PIj^V2 + PUvi) +■■■] 
« 1 + ({/*) ' [v^ - 2$Ar«2 + 2 {PIj^V2 + PU^i) +...] . 



* In the first post— Newtonian approximation tlie orders of tfie metric coefficients tliat are needed are 0(c~*), 0(c~^) and 
0(c~^) for goOi 90a and g^p respectively. However, as reported in Table 1 in CN, to determine the total energy without any 
additional assumptions we need g^^p to 0(c~*). An alternative method consists in evaluating the spatial metric gij to 0(c~^): 
in this case it is still possible to derive the conserved energy, provided that the equations are supplemented by a condition for 
iscntropic flow (see Ch65a, Chandrasekhar 1969b, c for a detailed discussion about this point). 
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We can proceed via successive approximations, substituting in the right hand of (|9|) the expansion 
(f7*)^ w 1 + (?7*)^ + . . ., where (L/*)^ is the (unknown) Newtonian contribution to (J7*)^. This yields 



(aUf 



1 + ft/* 



N 



[v^ - 2^NV^ + 2 (/3|,jv«2 + /3pjv«i) + 

1?^ 



9 / 



(92) 



At Newtonian order we have 



(93) 



which gives {U^)% — — 2$7v, and 



M3\ 



where (^ViP^) = /Jp^wi + /3p^W2- We obtain 



(94) 



The two expansions (|9^) and ( ^5| ) can be used together to yield 

7 



V 



V 15 o 7 n 3 ^ 

1 + y - 3$jv - 3$pjv + y - 2*w« + 8« 



and this expression can be written in a more compact form noting that 

5 



<^PNd^x = 



V 



V 



N 



d^3 



(95) 



(96) 



(97) 



The latter result can be verified by substituting equations (|73|), ([76|) directly and integrating; the answer 
has been checked by making use of an algebraic manipulator (MAPLE). As a consequence, an equivalent 
expression for the baryon rest-mass is 



45. 



13 



(98) 



The second conserved quantity that enters in our calculation is the total mass-energy M, which we 
take as the ADM mass (see e.g. Bowen & York, 1980) 

M = — - i V'-^d^S, = -— I \lHd^x. (99) 



27r J a 27r 

In equation (|9^) the surface integral is over the sphere at infinity, while Gauss's law has been used in the 
second step. Using the Hamiltonian constraint (pq), the mass can be rewritten as 



M = h+h, 



(100) 
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where 



To evaluate (101) we note 



while, by making use of expressions (BS), (r3) we find 



P~ Po 



Po 



o 



Af3 



(101) 



(102) 



(103) 



O 



(104) 



where y^i = dy/dxi. For simplicity, the subscript "PN" has been dropped and we simply write /3' for the 



leading term /3p^. The explicit form of (104) in terms of the index symbols is tedious but straightforward 
and is obtained by differentiating the shift function given by (p7|). The resulting expression has been 
computed with MAPLE. By collecting these quantities, exploiting ( |97| ) again and dropping from the 
integrand all terms of higher order, we obtain 



M 



Po 



2 4 Po 



PN 



1 



167r jv 

For later applications it is also useful to write down the conserved energy E = M — Mq 



(105) 



M - M, 







Po 



16tt 



PN 



(106) 



Finally, we need to evaluate the total angular momentum, J. This quantity is obtained from the 
integral (see e.g. Bowen & York, 1980) 



1 

16^' 



(107) 



where e^fe is the Levi-Civita tensor and integration is over the sphere at infinity. Let us also introduce the 
symmetric tensor K^^ = ^^K'^K Since it is 4" ^ 1 + O as r ^ oo, K"^^ may be used in place of the 

extrinsic curvature in (107), to lowest order. Moreover, again following Bowen & York (1980), this tensor 
satisfies 



kl 



kl 



(108) 
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It follows that 



IQtt 



(109) 



and, using (|108| ) and the momentum constraint (^7|), this can be rewritten as 

J, = e,,, j x^S^-^^H^x. 



(110) 



For rotation about the cca axis, we have J3 = J and the total angular momentum turns out to be 

J = csjkj x^S''^^°d^x (111) 



Exploiting equations (|l|), (|5|), (|6|) and (H), we are able to derive the PN approximation of the integrand 
in dml): 



5* 



po ( 1 + — ) (1 + 77" 



1 + 2« " 



{PI 



PN 



P 



pov' ( 1 - $Ar + 1'^ + — ) + Po/3pAr , 



and 



AT 



1 - G^N + V + 



(112) 



(113) 



Substituting (113) into (111), we find 



Po 



PN 



PN 



(Tx. 



(114) 



where all terms of order O [M^ / R^^ have been neglected in the integrand. The kinetic energy is T = VlJ /2 
[see CST92, equation (56)], and it is easy to recognize that (114) is correct at the PN level. 

The derivation of conservation laws in general relativity has been considered by Chandrasekhar 
(1969b) and CN, who used the symmetric energy-momentum complex of Landau-Lisfhitz to determine 
integral forms of conserved quantities at various post-Newtonian orders. Although CN do not present 
any expressions for the integrated conserved quantities, their results provide a useful tool to check the 



correctness of ( |106| ) and (114) by comparing the corresponding integrands. This comparison is presented in 
Appendix A. 3. 
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4-3.2. Evaluating The Conserved Integrals 
Summarizing, integral expressions for the conserved quantities that we have derived are 



V 



5 35 P 



PN 



Mo 
M- Mo 



1 

16^ 

Po 

Po 



K,jK'^d^x, 



15 



d^x, 



PN 



d^x 



(115) 



In order to perform the quadratures over the fluid volume, we adopt the ellipsoidal approximation for the 
matter distribution, whereas po — constant inside a triaxial ellipsoid. With this assumption, it is convenient 
to introduce and evaluate the following integrals: 



2:1 



Pod x = Mc = —^0010203 = '^Po^ ' 



6Ml 
5 R 



12 = po^Nd^x = 2M = 2Y^M^^ 

Jv 

13 = [ PQV^d^x = n^ [ po(xl+xl)d'^x = n^{In+l22)-^lMcn^R^^ 

Jv Jv 5/1 



(116) 
(117) 
(118) 



where A4 is the Newtonian potential energy, while R and the two dimensionless ratios / and h have been 
introduced in section (3) [see equations (Q), (0), In deriving the previous expressions we made use of 

the results 



i 

Mti = -2Trpo All u , 
Li. = 



PoXiXjd X = -McttiSij , 



(119) 

(120) 
(121) 



given by expression (2.4), (2.12), (2.13), (3.128), (3.129) in Ch69. Note that, in contrast to the Newtonian 
analysis, Mc is only a coordinate quantity and has no physical meaning. The other integrals appearing in 
equations ( 115 ) are 



Is 
2^10 



Po'^PN^'x , 15^ I Po'^^d^x , 

V 



(122) 



Pv^d^x, It ^ I po^NV^d^x 



..2j3„ 



/ pov'^d^x, 19= po^{xif3pj^ - (3pj^X2) d^x , 
Jv Jv 



'V 

1 

16^ 



K'^K.jd^x. 
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From the quadratic forms of $7v, ^pn, etc... it is easy to recognize that the evaluation of X^, i — 4. .9, only 
involves quadratures of the kind lij and Eij , where^ 

E^J = [ Pfix1x](fx = ^Mcala] (1 + 2,5^) ; (123) 
Jv •jo 

the same finding holds also in the case of /lo, but the integration involves some lengthy algebraic 
calculations. After many simplifications, and exploiting the properties of the index symbols, we obtain with 
the help of MAPLE the following integrated forms of the conserved quantities 

M « M, + 3ii^/+^M,02i?4 + i||.gi + il^f]2^V, (124) 
no a K'^ K 

Mo « Kh+^^^f + Ul^n^R^l + ^g. + ^^l'R^p,, (125) 
5 it 5 h R 

J « nM,R'-^(^l + ^^p3h^ , (126) 

3 Af^ 1 1 Af2 

M-Mo « -|il^/ + iA417^i?2i + iM|g^, + il^f]2^V2, (127) 
5 R o h R'^ R 

where the functions gi, pi are defined as (see Appendix D) 

99 ,2 9 1 ^2 of /,„o^ 



1"2 



56 224AiA2^' 



1"2 



,4 



6/ , 57^ A1A2 , 69 Ai+^2 , 3 1 , ^2 ^ ^ 

= 77^+35^^T+70^^+70A^(^^-^^) ^^'^^ 

21/ 57 , A1A2 33 Ai +^2 3 1 , , , ,2 , ^ 

= 2oi+56^^^+56^^+140AlA^^^-^^) ' ^^^^^ 

27 / , 171 , A1A2 , 111 ^1+^2 , 3 1 2 
P12 = P^-P^--Y^Qh^^,A,^ + —^-^ + —^^[A^-A,) (133) 

+ 2^10 , 

6/ , 24^ A1A2 , I8 A1+A2 . 3 1 2 
= 57^+35^^^+35^^+140V^(^^"^^^ • 

Due to the complexity of the K^^ Kij contribution appearing into Af, we simplify the notation by leaving the 

f \ — 1 

term Xio = I —^^^R? 1 Xio, without substituting the explicit form in terms of the ellipsoidal variables. 



For an expression in terms of these variables, see Appendix C. The structure of expressions ( |124|) -( 127 ) is 
particularly convenient for performing the required variations, since the full dependence on the two axial 
ratios is contained in /, h and, for the PN contributions, in gi and pi. We have checked that our result 
agrees with the PN correction to the Newtonian energy obtained by ST for homogeneous spheres (Appendix 
B.2). 



^The integrals Eij are easily evaluated by introducing the variables x'^ = Xi/ai, such that the equation of the ellipsoid is 
Si (^0 ~ ^' ^'^'^ using spherical polar coordinates in the new system. 
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5. Equilibrium Configurations 



As we did in the Newtonian case (see section 3), we now construct sequences of axisymmetric 
equilibrium models. Each sequence can be parametrized by, e.g., the value M/Rg of the spherical, 
nonrotating member. These configurations yield a PN generalization of Maclaurin spheroids and were 
originally investigated by Chandrasekhar (1965b) by using the tensor virial method. As we will confirm, 
the effect of general relativity is to attribute to a star of a given eccentricity a larger value of Vt^ / (tt/jq)- 
Sequences of relativistic, numerical models at fixed A/q^/Po have been also published by BI (see also 
Bonazzola & Schneider 1974). 

We construct the equilibrium sequence by minimizing M or, equivalently, M — Mq, keeping fixed Mq 
and J. Due to the complexity of the expressions, part of the calculations here and in the next section have 
been done by using an algebraic manipulator (MAPLE). The procedure is as follows. First, we combine 
expressions (126) and (127) and rewrite the quantity M — Mq as a function of i?, Ai, A2, po and J. This 
gives 

-/ + 7T1-TT7^+ ^.912 + — 77T^V23, (135) 



M - Mo ~ 
where Mc = Mc (po, R), P123 = P12 



5 R ■' 4AI,R^ i?2-9i2+ ^ ^3. 
P3 and where we used the relation 



n^R-^ 



M^R^ 4 
J2 25h^ 



M'^R^ 4 

The equilibrium sequence is then determined by minimizing M 



2 R^^ 
1 - ^^Psh 



d (M - Mo) d {M - Mo) dR 



dR 



Mo according to 
-Mo) 



0, 



R— const 



(136) 



(137) 



where, here and in the following, the partial derivative with respect to Ai is taken holding constant J, po 
and \j , with i ^ j 

Note that, in contrast to the Newtonian case, the mean radius R is no longer constant. The variation of R 
can be obtained exploiting the constraint dM^ = 0, which gives 



Since 



dR 



Mn ~ Mc + 



expression (139) can be approximated as 



dR 



Mr 



9/ 



OMq 
dR 



18 Mc. 
5 R^-' 



R.— C07lSt 



dMo 
'dX^ 



5 .P 



R— const 



4 Mci?2 



5 .P 
12 M^R 



dh 
dX^ 



R.— const. 



(139) 



(140) 



(141) 



In axisymmetry, we set ai — 02 when performing the two variations (137). Since, in this case, the first 
derivatives with respect to the two axial ratios are equal, we need only consider the case i = I. To simplify 
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the notation, we introduce the foUowing symbols for the derivatives of a given function X 



R— const 



\dX1dX2J 



dX 
'8X2 

d^x 



R— const 



dxl 



R— const 



At the PN order, the variation (137) with respect to Ai is 



25 .P , 
4 McR^' 



dXi 5 R ■' 4 Mci?2 



i?2 



(1) 



■512 



(1) 



and this expression can be combined with (141) to yield 

25 J2 



5-' 



4 M3i? i? 

5 



(1) 



4 M3i? ^ 
15 J2 



(1) 



0, (142) 



(143) 



0. 



AM^R' 2 A/3i? •' 48 Af6i?2 

It is now more convenient reexpress the previous equation in terms of the gauge invariant parameter 
Vl^ / (ttpo). From equation (126) we obtain 



J2 



A/3i? npo 25/i2 



1 



By using this result, the equilibrium condition (143) becomes 
4h^f^^'> 20 M, 



3 hW R 



5 



(144) 



(145) 



{h?Pl23. 



,(1) 



3/i(i) 3 21/(1) 



4;j2 V" ' 4 /i"^^ ' 20 ' 20 /i 

where, in the right hand side, the first term represents the Newtonian result, while the remaining terms give 
the PN correction. Finally we set ai — 02, using expressions (D13), (D14). The derivatives of / and h are 
derived using expressions (A3), (A9) in LRS, while for the derivates of (712 and P123 we used results reported 
in Appendix D.2 (expressions D16| ). Doing this, we obtain the relation between the angular velocity and the 
eccentricity along the relativistic equilibrium sequence. By adopting the same notation as in Chandrasekhar 
(1965b), we write the resulting expression in the form 



ill 



(1) 



2Mr 



ai 



E{e) 



(146) 



where A = Ai = A2 = (l — 
parameter 2Afc/ai, 

E{e) = 



and where the strength of relativity is measured by the compaction 



20 



A3/2 



9'^ 



18 



3 /,2 
4A^('^P123)^^ 



(1) 

8 A 



(147) 



3 / 21/(1) 
40 A2 + 20~ 
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is a function only of the eccentricity of the spheroid, and 



= 2 Ml - 



(148) 



which formally coincides with the Newtonian expression (p5|). However, we emphasize that the eccentricity 
which enters in our PN formalism (and that of Chandrasekhar 1965b) is defined in terms of the ratio of the 
coordinate quantities 03 and Oi and is not a gauge invariant parameter. 

The PN sequences of equilibrium are reported in Figure |l| for different values of the parameter M/Rg, 
where Rs is the equatorial radius in Schwarzschild coordinates, and, in the spherical limit 



ai 



Ml M 




(149) 



[see Lightman et al. 1979, expression (4), page 422]. Note that this function has a maximum for 
M/Rs = (M/Rs)^^^ = 5/18, corresponding to 2Mjai = 3125/11664 « 0.268. It follows that our PN 
formalism can be used to investigate relativistic sequences up to a maximum value / Rs) j^ax ~ 5/18 ~ 0.28. 
As we can see from Figure |l|, we find that in the relativistic case the value of ri^/(7rpo) is larger than what 
Newtonian theory predicts for the same value of the eccentricity, confirming the Chandrasekhar's (1965b) 
results. Sequences of equilibrium can be described in an equivalent way by using, in place of Vl^ / [-rpq), the 
ratio 



T 

W\ 



{\VLJ + Mo-M) 



(150) 



which reduces to the ratio of rotational kinetic energy to gravitational potential energy in the Newtonian 
limit. Note that r/|W^| is gauge invariant for rigidly rotating objects, since M — Mo, J and are gauge 
invariant (they can all be measured by observers at large distance from the spheroid). By using the 
expressions of J, M — AIq and setting ai =02, we find, to PN order, 



T 

W\ 



npo 2 (2Ai + A3A3) 



1 + 



2Mc 5 
ai 2VA 



2 



gl2 

3./ 



1 

4/ TTPO 



1 



P12 



TP3 



(151) 



Xi — A2 — A 



Figure || shows the sequence of equilibrium corresponding to the maximum value [M/ Rs)„^^^ = 5/18 
together with the Newtonian Maclaurin sequence {M/Rs =0). As we can see, in the range we are 
considering, T/jW^j is insensitive to the parameter M/Rg. Figure || shows the relation between the two 
gauge invariant quantities Q'^/npo and T/|P^|, for the same equilibrium sequences as in Figure 0; squares 
mark the secular instability point and will be discussed in the next section. 

The function E{e) and the corresponding function Ech obtained by Chandrasekhar (1965b) are 
compared in Table and Figure ^. As it can be seen, the two expressions are in close agreement up 
at least to the value of e corresponding to the secular instability point in Newtonian theory, with a 
maximum fractional difference of a few percent. The discrepancy increases at larger eccentricity (30 % at 
the Newtonian dynamical instability point). This discrepancy may be ascribed mainly to our coordinate 
ellipsoidal approximation for the deformation and equilibrium shape. This assumption is only exact in the 
Newtonian limit. For e <C 1, the two functions can be expanded as 



E w 0.219e2 + 0.1156-^ + 0.059e^ 



E, 



Ch 



0.228e2 + 0.1146" + 0.042e6 



(152) 
(153) 
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We have compared our results with the fuUy rclativistic, numerical sequences presented in Figure 1 in 
the paper by BI. To compare those computations with our PN models, we have restricted our attention to 
the least rclativistic sequence in BI, which corresponds to IM^jax = 0.206 (7^ = 0.154 in BI). To make the 
comparison, we note that BI parametrize in terms of a proper eccentricity, which is defined as 



\x\—X'2—^ 



dx?. 



= 1-^ 



dx\ 



where dp and de are the proper radii in the polar and equatorial directions, respectively. Using 
@ for we find 



2M, 

Oi 

2M, 



15 





3 




















840^ 



1 



(154) 

(155) 
|) and 

(156) 



which provides the relation between the two definitions of eccentricity. 

Figure m shows the comparison between the BI values of il^/ (Trpo) and the two PN sequences derived 
by using E{e) and the Chandrasekhar (1965b) expression Ecin with 2Mc/ai — 0.206; Table ^ contains the 
values of / (vrpo), evaluated by using our expression, and the corresponding quantities by BI. In order to 
make a comparison between the rclativistic corrections, which is more useful, we define the correction as 
the difference between the rclativistic and the Newtonian value of VL^ / {irpo) [given by expression ( |148| )] and 
compare these corrections. Since BI does not present tabulated values for Vl^ / (ttpq), the values we have 
reported in table ^ have been read off Figure 1 in BI, with an estimated error of ~ 0.02. For e < 0.7, this 
error is of the same order as the PN correction, making the measured value of the latter questionable. For 
this reason we have reported in Table || only the values of this quantity for e higher than 0.7, together with 
the corresponding fractional error. As it can be seen, in this case our values are lower with respect to the 
numerical ones by BI, but the difference in the corrections is at most w 10% for e w 0.8. 



6. The Secular Instability Point 

The determination of the point of onset of secular instability proceeds as in the Newtonian case [cf. 



equations (26), (|28|)], whereby we must evaluate 

92 (M - Mo) 



92 (M - Mo) 



Ai — A2 — A 



d\2d\i 



Ai — A2 — A 



By using expression (137), the second derivatives can be written as 

(M ~ Mo) d ( d {M ~ Mo) OR \d (M - Mo) 



(M - Mq) 
8X20X1 



d \ d{M~ Mo) dR 
Wi \ dR 
d f d (M - Mo) dR 
dX^ { dR dXl 



dX, 
d{M- Mo) 



dX,. 



R— const 



R— const 



(157) 

(158) 
(159) 
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where the expression in brackets is given by equation ( |143| ), multipHed by M'^/R. After some algebra we 
find that, at PN order, 



[M - Mq) 



+ - 



Ml 
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5 J2 



5^ ^4M3i? 
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35 J 
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4 A/3i? 
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(160) 
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(161) 



Given axisymmetry, we can set oi = 02 when equating expressions ( |16C| ) and ( |161| ); the equivalence between 
the subscripts 1 and 2 will be made explicit at the end of the calculations. After this substitution is made 
all terms containing first derivatives cancel in the final relation; hence, for the sake of simplicity, they will 
be dropped in the following expressions. Imposing condition (157) yields 



-^V[h] = ^v[f]-^l-fV[f] + -^fV[h] + -^hV[f] 
AM^R 5 ^-^^ i? \ 5 ^•'^ 4Af|i?^ ^ ^ 2 M^R ^ 

25 J2 



(162) 
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hV [h] + V [512] + 



4 Mm 



where the difference operator T) is defined as the difference between the two second derivatives of a function 
with respect to the axial ratios 

2?[X] = x(i2)_x("^ 



Replacing J with according to (144), expression (162) becomes 



ill 



Setting now Ai = A2 = A gives 
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is a function of the eccentricity. The last step consists in evaluating this condition along the equilibrium 



sequence. Using equations (146) and (14S) we obtain 



2M, 



= 2Bii + 



C(e) , 



where 



Cie) 



P{e) - -E{e) 



(167) 



(168) 



and the leading term 2Bii formally coincides with the Newtonian expression (|3^). The equation for C(e) 
has been derived by using MAPLE, using again the derivatives given in LRS [expressions (A3), (A9)] and 
the equations reported in Appendix D. The value of the eccentricity at the secular instability point, e^ec, is 
calculated by equating the two expressions (146) and (167) and is reported in Table ||. Figure ^ shows the 
ratio r/|Ty| evaluated at e = esec^ as a function of the compactness parameter M/Rg. Squares in Figure |^ 
mark the instability as a function of the two gauge invariant ratios fl^/{TTpQ) and r/|W^|, and thus separate 
the regions of stable and unstable configurations. As we can see, effect of general relativity is to move the 
instability point to an eccentricity larger than what Newtonian theory predicts. This corresponds also to 
larger values of the two ratios ri^/(7rpo) and (see Figure ||), showing, therefore, that relativistic 

gravitation tends to stabilize a star against secular instability to a Jacobi-like, nonaxisymmetric bar mode. 



7. Discussion and Conclusions 

In PN gravitation the rotational velocity along the equilibrium sequence is not much different from 
that one derived in the Newtonian limit, although a star of a given eccentricity has a slightly larger value 
of f2^/(7rpo)- General relativity, however, is more crucial in influencing the stability properties. According 
to our treatment, the critical value of the eccentricity for the onset of the bar instability increases as the 
star becomes more relativistic, in the regime in which the PN approximation is valid. Both invariant ratios, 
ri^/(7rpo) and T/|M^|, also increase at the onset of instability above the values found in Newtonian theory. 
Gravitational radiation does not enter at PN order; it is present only at the 2.5 PN level and higher. The 
secular instability we have identified is driven by the presence of viscosity. That conclusion is consistent with 
our variational procedure, whereby the rotation was kept uniform and angular momentum was conserved. 
Viscous dissipation conserves angular momentum and drives a star to uniform rotation. Gravitational 
radiation dissipation conserves circulation (Miller 1974; LRS), not angular momentum, and does not 
maintain uniform rotation. The presence of a stabilizing effect due to general relativity on the Jacobi-like 
bar mode instability is in agreement with the BFG and Bonazzola, Frieben & Gourgoulhon (1997) analysis 
of relativistic polytropes, who observe a growth of the critical adiabatic index when the relativistic character 
of the configuration is increased. Shear viscosity can provide a such dissipation mechanism in a very cold 
NS (T ^ 10^ K), but it is inefhcient for hot, newborn objects with T ^ 10^" K (see e.g. BFG), where the 
nonaxisymmetric instability is more probably induced by gravitational radiation (the CFS instability) and 
proceeds via a Dedekind-like mode. In this case, the fully relativistic SF computations show that effects of 
general relativity are reversed, and that the instability is significantly strengthed with respect to Newtonian 
theory. SF calculations are fully relativistic, so this bimodal behaviour may be explained in terms of strong 
field effects. However, as originally suggested by SF, there is also the likelihood that, in general relativity, 
the viscosity driven and the gravitational radiation driven m — 2 modes may no longer coincide and that 
the gravitational field plays a different role in each of them. Although at present no firm conclusion can 
be reached, this speculation is further strengthed by the results of our PN treatment. Unfortunately, a 
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direct comparison between our results and the ones obtained by BFG and SF is not possible, since in both 
that cases numerical constraints prevented them from treating incompressible fluid configurations, such as 
the Maclaurin and Jacobi ellipsoids, where the density profile is strongly discontinuous at the surface. In 
the BFG's case, numerical tests showed that the density steepens dramatically (and the numerical error 
increases more and more) for 7 > 3.25 {n < 0.44). Models presented by SF are restricted to values of 
the polytropic index n > 1, since the presence of discontinuous derivatives of energy density and metric 
functions across the surface of the configuration make the description less accurate for stiffer equations of 
state. Moreover, the m = 2 mode is only present in the SF model with n = 1. 

In Newtonian theory, the variational principle also has been used by LRS to investigate the secular 
instability to dissipation which conserves circulation C, rather than J, such as the emission of gravitational 
waves. In this case they used the energy functional for Dedekind ellipsoid, rather than a Jacobi ellipsoid 
as in the case of a viscosity-driven mode. Because the energy function of a Riemann-S ellipsoid is 
symmetric under interchange of C and J, the two analysis are virtually identical in Newtonian theory, with 
C appearing in place of J in all results. The previous discussion suggest that this symmetry may be broken 
in a relativistic treatment; this may be apparent at PN order and we plan to investigate this issue in a 
forthcoming work. 

We are grateful to J.L. Friedman and N. Stergioulas for providing a copy of Stergioulas's Ph.D Thesis 
prior to publication. We thank also Thomas Baumgarte for several very helpful discussions. This work was 
supported by NSF Grants AST 96-18524 and NASA Grant NAG 5-3420 at the University of Illinois at 
Urbana-Champaign. 



A. Comparison with the Chandrasekhar's Results 

Results presented here have been derived by using a 3+1 decomposition of the Einstein field equations. 
It is therefore useful to compare the resulting expressions with those obtained by CN, who worked with the 
usual field equations. In the first PN approximation, the metric coefficients goo, goa and gap are retained 
up to the order 0{c^*), 0{c^^) and 0(c~^), respectively (see CN and footnote ^ in this paper). In this 
Appendix we report the comparison. In section A.l we show that the same gauge condition is common to 
both treatments. In section A. 2 we compare the metric functions, while in A. 3 we compare the conserved 
quantities. 



A.l. The Spatial Gauge Condition 

Here we verify that the gauge condition adopted in CN is identical to our own. Chandrasekhar's 
expressions are obtained in the gauge 



where Pi — —Pppj [see CN, equations (3) and (6); section A. 2. 2 in this paper]. In our approach the 
condition is [see CST96, equation (3)] 

- Dkf3' , (A2) 



- 29 - 



where 7 = At leading order, we have the approximate expressions 



dxk 



dt 



7 dt 



6 



dt 



(A3) 



Substituting ^' 1 — $Ar/2 and neglecting terms of higher order, we obtain 

.9$ AT 



- 3- 



dt 



, dxk ' 



(A4) 



which coincides with the Chandrasckhar's choice (Al) 



A.2. The Metric 

A. 2.1. The Lapse Function and the Conformal Factor 

We can now compare expressions ( |63|) and (joj) for and apjy with the corresponding terms derived 
by CN. We start from the CN's expression of goo 

500 « 1 - 2C/ + 2 - 2$c/.) , (A5) 

where ^ch is the solution of 

VHch = -47rpo + [/ + ^— ) (A6) 

and t/ = — $Ar [see CN, equations (3) and (4)]. Note that CN used a different signature for the metric, i.e. 
H . At our order of approximation, we find 



ach ~ V9^~l-U + U^-2^ch-^ (A7) 
« 1-C/+ — -2$ch. 



On the other hand, ^pn and QpN are defined as the solutions of [equation (58), (|6l|)] 

V^^PN = -lOirpo^N + "^T^Pav^ , (A8) 

V^GpAT = 67rpo ( + 2— - i$jv ) , (A9) 



po 2 

so that, by comparing the corresponding differential equations, we can immediately recognize that 
—^Ch = ^PAr/4 + QpnI'2- Substituting this result into (A7) gives 

aCh « 1 + $Ar + ^ + ^ + QPN , (AlO) 

which agrees with the solution derived in section 4.1. 

For the spatial metric the comparison is trivial. The PN expression in CN is [CN, equation (3)] 

5afc « - (1 + 2C/) , (All) 

and, and the PN order, = -gab « (1 + U/2) 5ab- 
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A. 2. 2. The Shift Function 
Consider now the shift vector. In the CN's formalism, at the PN order it is [CN, equations (3), (5)] 

90a « « Pa , (A12) 

where Pa is the solution of the equation 

V^P.^-16wa + 0^. (A13) 



Using the gauge condition [CN, equation (6)] 



dU I ^ dP, , , , 



and neglecting terms of higher order, expression ( [A13| ) can be rewritten as 



3 dxa ^ dxi 

On the other hand, in our formalism we have [equations (^), (^4[), ([70|), (|7 



V2/3' = V^G" - (V'B) (A16) 



Taking the divergence of 13'^ and using equation (H) again gives 



4 4 

which can be substituted into (A_16), yielding Pi = — /3p^. 



= VfcG^- - iv^B = ^VfeG^- , (A17) 



A. 3. T/ie Conserved Quantities 



The derivation of conservation laws in general relativity has been considered by Chandrasekhar (1969b) 
and CN, who used the symmetric energy-momentum complex of Landau-Lisfhitz to determine integral 
forms of conserved quantities at various post-Newtonian orders.H In this Appendix we compare our results 



(106) and (114) with the corresponding integrands obtained by CN. Consider the conserved energy. In the 
first PN approximation, CN derive the energy per unit volume of the fluid [see CN, equation (67)] 



E 



Pa 



Po 



'-u^-lu- 
2 2 



n 



(A18) 



2 2 



2UII 







(-£) 









^Actually, in the first PN approximation the same expressions were originally obtained by Ch65a from a direct inspection 
of the PN generalization of the equations of motion and assuming the supplementary condition for isentropic flow. However, 
the fact that it is possible to derive the conserved quantities without the aid of the Landau-Lifshitz complex (at least when 
the additional condition is provided) is a peculiarity of the first PN order. The method discussed in Chandrasekhar (1969b) 
appears to be the most convenient way, allowing the derivation of the integral quantities in all the PN approximations beyond 
the first. 
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where 11 is the internal energy and Pi = — {Pi)pj~^- In the incompressible case of interest here, 11 = 0. 
In order to compare (A18) with the integrand function appearing in our expression (106), we note that 
CN's results have been derived exploiting a set of identities between quantities that are equal, modulo 
a divergence. In fact, upon integration, two functions that are equal modulo divergence give the same 
conserved quantities (see Chandrasekhar 1969b for a detailed discussion). After some lengthy reductions, it 
is possible to demonstrate that the K^^ Kij contribution appearing in ( |106| ) can be rewritten as 



1 / dPadPa ^ dPadPb 

2 \ dxb dxb dxh dxa 



From equations (62), (63) of CN we have 

dPa dPa 



dxb dxb 

dPg dPb 

dxb dxa 



IGirpoViPi - 3 



dU 
'dt 



i(divP)^ 



(mod div) . 



9 



dU_ 
'dt 



(mod div) , 



and these equalities can be used, together with the gauge condition (Al) 



to rearrange expression (A19). This yields 

K'^K,, EE 8ttpov,P' 



1 / -\2 1 / 

- [divPj - - [divPj (mod div) 



(A19) 

(A20) 
(A21) 

(A22) 

(A23) 



Sirpo [viP^)pj^ (mod div) 



Finally, inserting back (A23) into (106) one recognizes that our expression for the energy density agrees 
with equation ( A18 ). 

Consider next the angular momentum. Start from expression (141) by Ch65a for the conservation of 
total angular momentum, 

{xiTTj — XjiTi) (fix = constant , (A24) 



where 



TTi = avi + ipo {Ut - Uk-ik) + 4po [viU - Ui) , 



Ui, Uk-ik are quantities which enter the shift function, a = po (l 



2U + P/po), and where 



incompressibility implies 11 = 0. Specializing to our velocity profile, we have 7r3 = and 



Jch 



{xiTT2 - X27ri) tfx . 



(A25) 



(A26) 



After some algebra it is possible to demonstrate that the PN expression of the shift function in Ch65a is 



and that the integrand appearing into ( A26| ) can be rearranged as 
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4C/. + - (U, - Uk-ak) 
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rr ^ TT TT ^ rr 

-P0X1U2 - -poXiUk-2k + -j^P0X2Ul + -poX2Uk-lk + 



(A27) 
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Rewriting the last row in the previous formula as 



pi^xif? - poX2l3i = ("2/3^ + 



(A29) 



and inserting (A28), (A2£) into the conserved quantity (A26), we finally recover the same expression as in 



B. The Spherical Limit 

B.l. Comparison with the Exact Solution 

In the spherical limit, the interior metric for a homogeneous configuration admits an exact solution 
in full general relativity. In this Appendix we verify that our results agree with the PN expansion of the 
exact solution. We start from the interior 3~metric, written in Schwarzschild and in conformal (isotropic) 
coordinates as 



1 



"'^^^ - l-2m(r,)/r. 



where m{rs) = -iirpor^/S. Equating the metric coefficients gives 



1 



Vl-2m (r,)/r 



zdr. 



and these relations can be combined to yield the differential equation 



-d* = 1 



1 



\/l-2m {rs)/rj r. 



Integrating (B4) gives 



* = if 1 



1 



1/2 



(Bl) 

(B2) 
(B3) 

(B4) 
(B5) 



where K is a, constant of integration. By using (B3), expression (B5) can be written in terms of isotropic 
coordinates as 

(B6) 



* = 



1 + 



1/2 



The value of K is determined by matching the solution with the exterior metric. This is [see Misner, Thorne 
& Wheeler, equation (31.22), page 840] 

M 

vl. = l + -, (B7) 
and, at the surface of the star, r = R. Imposing the boimdary condition yields 



V2K = { 1 



M 
2R 



3/2 



The PN expression of ^ can be then derived expanding (B6), and using the relations 



(B8) 



(B9) 
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This yields 



+ (BIO) 

* « 1 + npoR' - + f {npof ~ \ {npof + ^ {npof . (BU) 

3 6 3 6 

By specializing our solution ( [73|) and (|7^) to the spherical case gives 

$ = + a>p^ = -27rpoR' + - | {npof R' + f {npof R'r' - ^ {npof • (B12) 



which, with (59), agrees with (Bll). 



B.2. The Relativistic Correction to the Total Energy 

Our integrated expression for M — AIq allows a comparison with the relativistic correction to the total 
energy of a spherical configuration as derived by ST (see also Zel'dovich & Novikov, 1971). Following ST, 
let us define the relativistic correction as 

AEgtr = M - Mo - En , (B13) 

where En is the Newtonian gravitational energy 



5 Rn 5 V 3 

and Mjv, Rn are the Newtonian mass and radius, respectively. Dimensionally, at PN order this correction 
can be written as 

AEgtr = -kM^^^pT . (B15) 

and the value of k in the homogeneous case can be obtained by using equation (6.9.29) in ST for an 
incompressible gas with polytropic index n = 0. This yields]] 

3 /47r^^/^ 



In order to recover this result from our equations, we specialize equation (127) to the nonrotating, spherical 
case. In this limit Ai = A2 = = 2/3 and / = 1, so that 



'Note that, in the Hmit n — > 0, the two integrals appearing in ST (6.9.29) become h = -3M^/{14:R%) and h 
-3M3/(35i?|f). This gives 

3 A/3 3 /47r\2/3 2/3 
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On the other hand, the exact relation between the conformal (isotropic) radial coordinate R and the 

(B18) 



Schwarzschild radius Rg is given by ( BIO ). At PN order 



R^RAl- 



M 



Rs 

and we left to derive the relation between Rg and Rn- Start from the definition (6.9. f2) in ST 



1/3 



where V is the proper volume. To first order [see (6.9.f 5)] 



Rn 



Rs 1 



Rs 1 



^s JO 



mrgdrg 



(B19) 



(B20) 



This gives 



1 47rpo 
Rs 3 



Rs ~ Rn I 1 



M 
5R^ 



Rs 1 



M 
5Rs 



and, combining (Bf 8) and (B21), yields 



R ~ Rn 1 — T 



6 M 



5R 



N 



(B2f) 



(B22) 



Finally, substituting the equation (B22) into equation (B17) and subtracting En, we obtain the relativistic 
correction 

3 Af3 



M - Mo ~ En 



3 Af|f 
'70 Rjf 



70 R% ' 



(B23) 



which agrees with the ST's result. 



C. The Integral Iw 

The integral Xio has been obtained by using MAPLE. We only report here the resulting expression, 
that is 

iio = (39/560)(Ai-^2)V(AiA2)-(9/f40):4i2(f-A?/A3)(Ai22-Aii2A3/A3)A2/A^ (Cf) 

+ (3/2240)(Ai - ^2)2(9^112 + AI22 + {-'2A2XI + 2X1 - A-iXl 

- A12 A? + A12 A?A^) V(f - Xir/Xt/{1 - Xlf)X\/Xl 

+ (3/8)(A?/A4)(-f/5(^i - A^XllXlf + (3/35)A?(-f /(f - A?)(2Ai - 2 + A2) 

+ f/(f - Xl){2A2 - 2 + A^)f + (6/35)A3(-f /(f - A?)(2Ai - 2 + A2) + f /(f - Xl) 

X (2^2 - 2 + Ai))(Ai - (f/3)3i2(A3/A^) + f/3/(f - Xl){2Ai - 2 + A2)Xl 

- {Xl/Xl){A2 - (1/3)^12 + f/3/(f - Xl){2A2 - 2 + Ai)Xl))) + (3/280)^i2(f 

- A?/A3)(A?/A^)(33i2(Ai - A2XI/XI) - 9^1^112 

+ %A2{Xl/Xl)Ai22 + 3(Aii2 - Al22X\/Xl)Ai2Xl/Xl 
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- A?/(l - Xl){2Ai -2 + ^2)(3Aii2 - i-2A2\l, + 2X1 - ^i^i - 3i2A? 

+ Ai2XlXl)/i-l + Xl)/{-l + XD) + Xl/{1 - Xl)i2A2 - 2 + ^i)(33i22A?/A^ 

- (-2A2A3 + 2A3 - AiAf - Ai2Xl + A,2XlXl)/i-l + Xl)/{-l + A?))) 

+ (3/70)(A^/Ai)(l/(l - A3)2(2^2 - 2 + A^ fxl/Xl + 1/(1 - Xff{2Ai - 2 + ^2)^) 

+ (3/1120)^12(1 - A?/A^)2(l + Xl/Xl){-2A2Xl + 2X1 - ^1^2 - 3i2A? 

+ A,2XlXlf/{-l + A3)VAi/(-1 + A?)VA^ - (3/140)^12(1 

- A?/Ai)(-1/(1 - A?)(2Ai - 2 + A2) + A?/A3/(1 - A3)(2A2 - 2 + Ai)){-2A2 * A^ 
+ 2Ai - AiXl - A12XI + Ai2A?A2')/(-l + A^)/Ai/(-l + A?)/A2 

+ (3/128)3i2(l - Xl/Xlf{Xl/Xi){-Al2 + (6/7)Ai2An2 + (9/35)3ii2 

+ (2/7)Ai2(-2^2Ai + 2Ai- AiAi-3i2A? 

+ 3i2A?A^)/(-l + Ai)/(-l + A?) + (l/35)(-2A2A^ + 2A^ - A^Xl - 3i2A? 

+ 3i2A?Ai)V(-l + Ai)V(-l + A?)' + (27/35)]4i22A?/Ai - (6/35) A122 (-2^2 A^ 

+ 2A2^ - A,Xl - A12XI + 3i2A?Ai)/(-l + A^)A?/(-l 

+ A?)/Ai) + (3/10)Ai2A?/A^ + (81/280)(A2/At)(Ai - (1/3)3i2A?/A3 + 1/3/(1 - A?)(2Ai 

- 2 + A2)Xlf + {81/280){Xl/Xi)iA2 - (l/3)Ai2 + 1/3/(1 - A3)(2A2 - 2 + Ai)Xl)^ 

- (27/140)(^i - (1/3)3i2A?/A3 + 1/3/(1 - A?)(2Ai - 2 + A2)Xl){A2 - (l/3)Ai2 + 1/(3 

- 3A3)(2A2 - 2 + Ai)A3)/Ai/A2 - (9/70)3i2(Ai + A2A?/A3)/Ai/A2 + (3/70)Ai2A? 
X (-1/(1 - A?)(2Ai - 2 + A2) - 1/(1 - Xl){2A2 - 2 + Ai))/X2, 

where 

A12 = af^i2 , A122 = 0^^122 , A112 = 0^^112 , (C2) 
only depend on the axial ratios. 



D. The Index Symbols 

The expressions of the metric functions we have derived have been presented in terms of the index 
symbols Aijk... and Bijk..., where 

Jo A{al+u){a] + u){al + u)... 

Jo A{al+u){a^j+u){al + u)... 

and where = {a1 + u){a2 + u){a^ + u). The use of these quantities is particularly convenient, since it 
allows for a more compact form of the results. In the expressions of the conserved quantities, however, it 
is more convenient to reduce the number of different index symbols entering in the final forms. This is 
because when we perform the variations of the conserved quantities by using MAPLE, the derivatives of 
the index symbols are evaluated analytically and then substituted in the variations. The index symbols 
are not independent: they are manifestly symmetric in their indices and obey to a set of identities and 
recursive relations (see Chapter 3 in Ch69). Exploiting these properties we derived the results presented in 
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D.l; these have been used to obtain the final expressions of the gi and pi functions which enter in A/, Mq 
and J [see (128)-(134), (CI)]. Note that the final forms only depends on the quantities 

A,,Ai2,Ai22,Aii2 . (D3) 

Then, variations have been performed by using a MAPLE program. The final results have been specialized 
to axisymmetry, using expressions (D13), (D14) and substituting the derivatives given in D.2. 



D.l. The Coefficients Aijk and Bijk; General Expressions 

The Ai's functions are given in terms of the standard (elliptic) integrals by expressions (3.33)"(3.35) in 
Ch69. The other index symbols can be written as 

Aa = — j^{A,-Ai) iovi^l (D4) 

" ~ Zaf 3 3 



2 -^a " o Aim ioT i^l^m 



3 

1 



aj ~ at 



{A,m - Aim) fOTi^l (D6) 



A^ll = (^«' " ^u) iori^l (D7) 



a: - at 



oa^ af 3 a. 

5 

3a. 



2 Ail 2 o T for i ^ Z 7^ m 

Ay^^»'™«m(l-'5mi)(l + 2<5/m) toi i I (D8) 

An. ^ — ^ 



1 af o? 

1 

5^ 



The i?i;m..'s are then: 

B,=I~ ajA, , (DIO) 

SiJm.. = Aim.. — afAiim.. ■ (DH) 

By using these relations we obtain the final forms of pi , qi and Xio , which only contain the quantities 

A,Ai2,Ai22,An2. (D12) 

Moreover, the A^'s are not independent, since A3 = 2 — Ai — A2. The derivatives of these quantities have 
been evaluated analytically and are reported in the next section. In axisymmetry, we have A3 = 2 — 2Ai 
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and 



Al2 
Aii2 



Au 



5e2 



4e4 



(3Ai - 2) 



(D13) 
(D14) 



D.2. Some Useful Derivatives 

The first and second derivatives of tlie quantities At, A12, A122, ^112 with respect to the axial ratios 
have been evaluated using the chain rule 



dA. 



ijl 



E 



dAiji,, dak 
dak dX^ 



(D15) 



and exploiting the definitions (Dl). The resulting expressions have been evaluated in the axisymmetric 
case, and substituted into the variations of the conserved quantities by using MAPLE. The final results, 
already specialized to axisymmetry, are 



and 



3X2 
dXi 

dAi2 
8X2 

dAi22 

dXi 

dAi22 

8X2 

dAn2 
dXi 

dAii2 
dX. 



dAi 


15 Ai + l2AiX^ - 


18 A3 


dXi 


(8 -8 A3) A 




dAi 


-3Ai + 12yliA3 - 


-6 A3 


8X2 


(8 -8 A3) A 




dA2 


-3Ai + 12AiA3 - 


-6 A3 


dXi 


(8 -8 A3) A 




dA2 


15 Ai + 12 AiA3 - 


18 A3 



(8 -8 A3) A 
9 + 18 A3 + 12 A*^ - 54 ^1 A3 
16 {\-X^fX 

27 - 42 A3 + 12 A^ + 18 A3 

16 (1- A3)^A 
-615 Ai + 306 A3 - 1336 A^ + 400 A^ + 1560 AiA3 
384 (1 - A3)3 A 
-1095 Ai + 2034 A3 - 1448 A^ + 464 A^ - 480 AiA3 
384 (1 - A3)3a 
15 Ai - 738 A3 - 424 A^ + 112 A^ + 1560 v4iA3 
384 (1 - A3)3a 
-465 Ai + 990 A3 - 536 A^ + 176 A^ - 480 A3 



d^Ai 
dXl 



384 (1 - A3)3a 

39 Ai - 306 A3 + 516 AiA3 + 120 ^1 A^ - 144 A^ 
32 (1- A3)^A2 



(D16) 
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-9^1 + 72AiX^ - 18 A-"^ - 72 + 72 A1(A, 12)X^ 

32 (1-A3)^A2 
3 + 6 A3 - 96 A^ + 12 AiA^ + 120 AiX^ 

32 (1 - A3)2a2 
-9 -li + 72 .liA"' - 18 A-' - 72 A<' + 72 .11(A. r2)X^ 

32 (1- A3)2a2 

_ -117Ai + 144^1 A^ - 234A3 + 3048A6 + 336A9 - 4752AiA6 
256 (1 - X^ f A2 

(D17) 

-195 Ai - 960 AiX^ + 634 A^ + 1144 A^^ + 112 A^ - 1680 A^X^ 
256 (1 - X^f A2 

8325 Ai - 20220 AiX^ - 7926 A^ + 18624 A^ - 21296 A^ + 4928 A^^ + 20400 AiX^ 

512A2(1- A3)^ 

5535 Ai - 10440 AiX^ - 7362 A^ + 21528 A^ - 11184 A^ + 2688 Ai2 - 3600 AiX^ 

512A2(1- A3)^ 

-75y4i -480y4iA3 + 1386 A^ - 12424 A'' - 2704 A^ + 20400 AiA^ + 512 A^^ 

512A2(1- A3f 

_ -345 Ai - 4560 AiX^ + 2382 A^ + 3272 A^ + 16 A^ - 3600 A'^ 

512A2(1-A3)^ 
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Fig. 1. — The ratio fi^/ (ttpq) versus the eccentricity for PN equihbrium sequences of constant rest mass 
(soHd lines). The different curves correspond to 11 equally spaced values of the compaction parameter M/Rg 
in the range [0., 0.275]. This parameter characterizes the nonrotating spherical member of each sequence. 
The dashed line is the Newtonian Maclaurin sequence [M/Rg = 0). 

Fig. 2. — The ratio versus the eccentricity for the PN equilibrium sequence with {M/Rs)max ~ 5/18 

(solid line), and the Newtonian Maclaurin sequence (dashed line). 

Fig. 3. — The ratio O^/ {ttpq) versus T/IVFI for the same sequences as in Figure 1 (solid lines); the dashed 
line is the Newtonian Maclaurin sequence. Squares mark the secular instability point. . 

Fig. 4. — Comparison between the ellipsoidal PN function E{e) (solid line) and the corresponding function 
Ech derived by Chandrasekhar (1965b, dotted line). 

Fig. 5. — The ellipsoidal PN equilibrium sequence (solid line) and the Chandrasekhar (1965b) sequence 
(dotted line) for 2Mc/ai = 0.206. Squares represent the Butterworth & Ipser (1976) numerical values, and 
the dashed line is the Newtonian Maclaurin sequence. 

Fig. 6. — The critical ratio T/|VK| at the secular instability point versus the compaction parameter M/Rg 
along a PN equilibrium sequence (solid line). The dashed line is the Newtonian value. 



Table 1. The Function E{e); Comparison with Chandrasekhar (1965b) 



e E{e) Ech{e) err* e E{e) Ech{e) err* 














0.20 


0.0090 


0.0093 


0.0397 


0.25 


0.0142 


0.0147 


0.0389 


0.30 


0.0207 


0.0215 


0.0379 


0.35 


0.0287 


0.0298 


0.0366 


0.40 


0.0383 


0.0397 


0.0350 


0.45 


0.0480 


0.0513 


0.0661 


0.50 


0.0629 


0.0650 


0.0268 


0.55 


0.0790 


0.0808 


0.0220 


0.60 


0.0974 


0.0991 


0.0177 



0.65 


0.1191 


0.1202 


0.0094 


0.70 


0.1448 


0.1444 


0.0026 


0.75 


0.1752 


0.1719 


0.0195 


0.80 


0.2118 


0.2026 


0.0447 


0.8127* 


0.2223 


0.2108 


0.0531 


0.85 


0.2564 


0.2355 


0.0847 


0.90 


0.3115 


0.2666 


0.1552 


0.9529** 


0.3838 


0.2768 


0.3240 


0.96 


0.3939 


0.2723 


0.3651 



*The fractional difference is defined as 



\E{e)-EcH{e)[ 
- E{e) + Ech{e) " 



* Newtonian value for the secular instability point. 
**Newtonian value for the dynamical instability point. 
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Table 2. Comparison with Butterworth & Ipser (1976); 2Mc/ai = 0.206* 



es7 



ML 



ill 



BI 



ill 



AT 



nl 



AT. 



err 



BI 



0.199 


0.200 


0.024 


0.025 


0.002 






0.298 


0.300 


0.054 


0.058 


0.006 






0.397 


0.400 


0.097 


0.107 


0.011 






0.498 


0.500 


0.155 


0.164 


0.018 






0.598 


0.600 


0.227 


0.247 


0.027 






0.695 


0.700 


0.314 


0.329 


0.040 


0.051 


0.24 


0.796 


0.800 


0.419 


0.429 


0.059 


0.066 


0.10 


().84() 


0.<S.')0 


0.472 


0.474 


0.072 


0.0(i(i 


0.08 



'^Corresponding to 7s = 0.154, where 7s is the parameter used by Butterworth & 
Ipscr (1976). 

^The fractional difference with respect the Butterworth & Ipser value is defined 

as 

^Q~QBI\ 



err = 2- 



Q + Qbi 



where Q = ny {npo) - [OV (npo)] ^. 
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Table 3. Location of the PN secular instability point 



Rs ^««'= Wo J sec Vl^i; 





0.'^ 


0, 


.8127 


0, 


.3742 


0, 


.1375 


0, 


.025 


0, 


.8561 


0, 


.4230 


0, 


.1673 


0, 


.050 


0, 


.8850 


0, 


.4582 


0, 


.1924 


0, 


.075 


0, 


.9042 


0, 


,4831 


0, 


.2120 


0, 


.100 


0, 


.9171 


0, 


.5011 


0, 


.2268 


0, 


.125 


0, 


.9261 


0, 


.5144 


0, 


.2378 


0, 


.150 


0, 


.9325 


0, 


.5247 


0, 


.2458 


0, 


.175 


0, 


.9372 


0, 


.5327 


0, 


.2516 


0, 


.200 


0, 


.9403 


0, 


.5384 


0, 


.2554 


0, 


.225 


0, 


.9425 


0, 


.5424 


0. 


.2578 


0, 


.250 


0. 


.9436 


0, 


.5446 


0, 


.2591 


0, 


.275 


0, 


.9441 


0, 


.5457 


0, 


.2597 



^Newtonian limit. 



0.3 



0.28 - 



0.26 - 



0.24 - 



0.22 - 



0.2 - 



0.18 



0.16 



0.14 








0.05 



_M_ 



0.1 



0.1 



